Summary work

Setting up the environment

# Packages
library(ggplot2)
library(dplyr)
library(plotly)
library(tidyr)
library(RColorBrewer)
library(limma)
library(gplots)
library(edgeR)
library(devtools)
library(ggbiplot)
library(DESeq)
library(geneplotter)
library(Rtsne)
set.seed(1234) # makes everything reproductible 
meta_data <- read.table("../Data/raw_data/meta_data.txt", header = TRUE)
dataInitial <- read.table("../Data/raw_data/expn_matrix.txt", header = TRUE, row.names = 1)
spike_data <- read.table("../Data/raw_data/expn_matrix_spike.txt", header = TRUE, row.names = 1)

Normalization

Kevin’s normalization

Here the data is normalized based on the spike_data we have. Although this seems to be correct, there could be some issues: See post.

spike_data <- rbind(spike_data, "Totals" = colSums(spike_data)) #finding sum of column and adding it as new row
scaled_totals <- spike_data["Totals", ]
average_reads <- rowMeans(scaled_totals)
scaled_totals <- rbind(scaled_totals, "Scale" = scaled_totals[1, ] / average_reads) #dividing total spike values by average spike values
dim_names <- list(row.names(dataInitial))       #setting dim_names for upcoming matrix
dim_names[2] <- list(colnames((dataInitial)))
normalized_data <- matrix(nrow = nrow(dataInitial), ncol = ncol(dataInitial), dimnames = dim_names) #pre-making a matrix
for(i in 1: ncol(dataInitial)){   #for loop which scales data and saves it as normalized_data
  normalized_data[ , i] <- dataInitial[ , i] / scaled_totals[2, i]
}
normalized_data <- as.data.frame(normalized_data)
#normalizing for cell size (cellular miRNA counts)
total_data <- rbind(spike_data, "Total Spike" = colSums(spike_data), "Total Reads" = colSums(dataInitial))
total_cellular <- total_data["Total Reads", ]
average_reads <- rowMeans(total_cellular)
cell_scales <- rbind(total_cellular, "Scale" = total_cellular[1, ] / average_reads) #dividing total spike values by average spike values
dim_names <- list(row.names(dataInitial))       #setting dim_names for upcoming matrix
dim_names[2] <- list(colnames((dataInitial)))
normalized_data2 <- matrix(nrow = nrow(dataInitial), ncol = ncol(dataInitial), dimnames = dim_names) #pre-making a matrix
for(i in 1: ncol(dataInitial)){   #for loop which scales data and saves it as normalized_data
  normalized_data2[ , i] <- normalized_data[ , i] / scaled_totals[2, i]
}
normalized_data2 <- as.data.frame(normalized_data2)
# deleting dead cells
alive_indexes <- list() #initializing lists
k <- 1  #initializing counter
for(i in 1:ncol(total_data)){
  if(total_data["Total Reads", i] > 20000){
    alive_indexes[k] <- i
    k <- k + 1
  }
}
alive_data <- normalized_data2[ , as.numeric(alive_indexes)]
alive_meta_data <- meta_data[as.numeric(alive_indexes), ]
#Saving completely normalized data as n_data and n_meta_data
n_data <- alive_data
n_meta_data <- alive_meta_data
#makes group and design matrix
data_day <- as.character(n_meta_data$Population)
data_day <- recode(data_day, "HL60D0" = "Day 0", "HL60D1" = "Day 1", "HL60D3" = "Day 3", "HL60D7" = "Day 7")
data_day <-factor(data_day)
design <- model.matrix(~0 + n_meta_data$Population)
#deletes "n_meta_data$Population" from name of columns
colnames(design) <- gsub("n_meta_data\\$Population", "", colnames(design))
#lets call intercept HL60D0
#colnames(design)[1] <- "HL60D0"
# Keeps genes without all zeros miRNA
dataNZ <- n_data[which(rowSums(n_data) > 0),] 
#making the contrast matrix
# let's suppose for now that we want to compare each group with each other ?????????????????????????
contrastMatrix <- makeContrasts(HL60D1-HL60D0,
                                HL60D3-HL60D0,
                                HL60D7-HL60D0,
                                HL60D1-HL60D3,
                                HL60D1-HL60D7,
                                HL60D3-HL60D7,
                                levels = c("HL60D0","HL60D1","HL60D3","HL60D7"))
#will not be needed: don't keep big data
remove(normalized_data)
remove(normalized_data2)
remove(alive_data)

Please note that with this normalization, we are only using cells ????????????????(alive cells?)

General normalization

To be sure that the normalization is good I will continue the analysis with a more general normalization. The following is based on this page. By general I mean not specific to pikes.

data_dayDESeq <- as.character(meta_data$Population)
data_dayDESeq <- recode(data_dayDESeq, "HL60D0" = "Day 0", "HL60D1" = "Day 1", "HL60D3" = "Day 3", "HL60D7" = "Day 7")
data_dayDESeq <-factor(data_dayDESeq)
deSeqDat <- newCountDataSet(dataInitial, data_dayDESeq)
# Note: actually it's not a real normalisation. Rather computing size factors depending on ratio of medians 
# if all size factors are roughly equal to one, the libraries have been sequenced equally deeply.
deSeqDat <- estimateSizeFactors(deSeqDat)
head(sizeFactors(deSeqDat))

HL60D0.C01_S1 HL60D0.C02_S2 HL60D0.C03_S3 HL60D0.C04_S4 HL60D0.C05_S5 HL60D0.C06_S6 
    0.5254052     1.1318130     0.9646220     0.3233751     1.3896389     0.9596427 

idx.nz <- apply(counts(deSeqDat), 1, function(x) { all(x > 0)})
nNZsamples <- sum(idx.nz)
#will not be needed: don't keep big data
remove(dataInitial)

We see that the number of non zero samples in all genes is very low compared to traditional RNA-seq: 5 nomrally in thousands. This may be a good reason to stick with the normalization with skie RNA.

#plotting the estimated dispersions against the mean normalized counts
deSeqDat <- estimateDispersions(deSeqDat)
plotDispEsts(deSeqDat)

multidensity( counts(deSeqDat, normalized = T),xlab="mean counts", xlim=c(0, 1000))

 multiecdf( counts(deSeqDat, normalized = T),xlab="mean counts", xlim=c(0, 1000))

The two charts above clearly shows that the second normalisation isn’t convincing. Indeed, the second chart assesses whether the normalization has worked, and the densities should overlapp since most of the genes are heavily affected by the experimental conditions. Note: The strange density chart could also be due to the fact that miRNA are very rarely expressed in every sample (as we have seen before).

Plots

PCA

PCA no standardization

First we use a log transformation to make the data approximatively follow the homoscedasticity assumption.

#Log transform
log_normalized_data <- log(dataNZ)
log_normalized_data[log_normalized_data == "-Inf"] <- 0
data_pca <- prcomp(t(log_normalized_data))
g <- ggbiplot(data_pca, 
              scale = 1, 
              obs.scale = 1, 
              varname.abbrev = FALSE,
              var.axes = FALSE,
              pc.biplot =TRUE,
              circle = TRUE, 
              groups = data_day, 
              ellipse= TRUE) +
  ggtitle("PCA without standardizing") +
  theme(text = element_text(size=14),plot.title = element_text(hjust = 0.5)) 
g 

#will not be needed: don't keep big data
#remove(g)
remove(data_pca)
remove(log_normalized_data)

Note: The PCA graph is interesting as it shows that after each day the clusters seem to go down in PC2 and extend their variance in PC1.

PCA standardized

Knowing the importance of having standardized data for PCA, we could think of standarding ours. This in general is a good idea, as veraibles are often not on the same scale and thus cannot be compared directly. Here, on revanche we have dimensions that are comparable as they each represent the number of miRNA in each cell. It is therefore a good idea not to standrdize, as it will unnecessarily decrease the variance explained. More information on this can be found on the following forums: (here)[http://stats.stackexchange.com/questions/105592/not-normalizing-data-before-pca-gives-better-explained-variance-ratio] or (here)[http://stats.stackexchange.com/questions/69157/why-do-we-need-to-normalize-data-before-analysis]. Out of curiosity we will still plot the standardized data:

#Log transform
log_normalized_dataNZ <- log(dataNZ)
log_normalized_dataNZ[log_normalized_dataNZ == "-Inf"] <- 0
data_pcaNZ <- prcomp(t(log_normalized_dataNZ), center = TRUE, scale. = TRUE)
g <- ggbiplot(data_pcaNZ, 
              scale = 1, 
              obs.scale = 1, 
              varname.abbrev = FALSE,
              var.axes = FALSE,
              pc.biplot =TRUE,
              circle = TRUE, 
              groups = data_day, 
              ellipse= TRUE) +
  ggtitle("PCA non zero genes") +
  theme(text = element_text(size=14),plot.title = element_text(hjust = 0.5)) 
g 

#will not be needed: don't keep big data
#remove(g)
remove(data_pcaNZ)
remove(log_normalized_dataNZ)

From this plot we see that as we’ve supposed, the variance explained is lower than in the first one. We thus conclude that we should stick with the first plot as we can (because on the same scale).

PCA top 50

Let’s now plot the top 50 genes

#arranging data based on total expression
indexTopFifty <- sort(rowSums(n_data), index=T, decreasing=TRUE)$ix[1:50]
topFifty <- n_data[indexTopFifty,]
# note: small difference with kenny probably come from the fact that he did log before and then took top fifty 
#Log transform
log_topFifty <- log(topFifty)
log_topFifty[log_topFifty == "-Inf"] <- 0
#PCA
data_topFifty <- prcomp(t(log_topFifty))
g <- ggbiplot(data_topFifty, 
              scale = 1, 
              obs.scale = 1, 
              varname.abbrev = FALSE,
              var.axes = FALSE,
              pc.biplot =TRUE,
              circle = TRUE, 
              groups = data_day, 
              ellipse= TRUE) +
  ggtitle("PCA top 50 genes") +
  theme(text = element_text(size=14),plot.title = element_text(hjust = 0.5)) 
g

#will not be needed: don't keep big data
#remove(g)
remove(data_topFifty)
remove(log_topFifty)
remove(topFifty)

Correlation heatmap

# Note: both plots of dataNZ and full data gives same correlation heatmap
normalized_data_cor <- cor(dataNZ)  #calculating sample-sample correlations
scaled_normalized_data_cor <- scale(normalized_data_cor)  #scaling for better contrast
limit <- 2
scaled_normalized_data_cor[which(scaled_normalized_data_cor > limit)] <- limit    #trimming for better contrast
scaled_normalized_data_cor[which(scaled_normalized_data_cor < -limit)] <- -limit
mypal <- rev(colorRampPalette(brewer.pal(11,"RdBu"))(100)) #setting colour palette
alphabet <- c("#b2df8a","#33a02c","#ff7f00","#fdbf6f","#b15928","#1f78b4","#6a3d9a","#e31a1c","#a6cee3","#cab2d6","#ffff99","#fb9a99")
pops <- alphabet[as.factor(data_day)] #setting population colour palette
heatmap.2(scaled_normalized_data_cor, Rowv = NA, Colv = NA, symm = T, #making a heatmap
          trace = "none", dendrogram = "none", col = mypal,
          labCol = NA, labRow = NA,
          cexCol = 0.45, cexRow = 0.45, main = "Sample-Sample Correlations",
          ColSideColors = pops) 

#will not be needed: don't keep big data
remove(normalized_data_cor)
remove(scaled_normalized_data_cor)
#remove(heatmapPlot)

In this heat map we clearly see the high correlation between the samples of day 0 and 1 and between the samples of day 3 and 7.

t-SNE plots

General t-SNE

Let’s plot the t-SNE plots of the cells that aren’t dead to see if we can find clusters and / or outliers. Note that we aren’t taking into account genes that aren’t expressed in any samples.

# t-SNE plot with alive cells and non zero genes
#ordering data by expn
alive_tsne_data <- cbind(dataNZ, "avg_expn" = rowMeans(dataNZ), "miRNA" = row.names(dataNZ)) 
alive_tsne_data <- dplyr::arrange(alive_tsne_data,desc(avg_expn))
alive_tsne_data <- t(alive_tsne_data)
colnames(alive_tsne_data) <- alive_tsne_data["miRNA", ]
#removes rows that were added for ordering
alive_tsne_data <- alive_tsne_data[!rownames(alive_tsne_data) %in% c("miRNA","avg_expn"), ]
tsne_out <- Rtsne(dist(alive_tsne_data)) # Run TSNE
# Show the objects in the 2D tsne representation
tSNETotal <- qplot(x=tsne_out$Y[,1], y=tsne_out$Y[,2],color=data_day) + 
  labs(colour = "Cell Type", x = "tsne1", y = "tsne2")+
  ggtitle("t-SNE alive cells and non zero genes") +
  theme(text = element_text(size=14),plot.title = element_text(hjust = 0.5)) 
tSNETotal

#will not be needed: don't keep big data
remove(tsne_out)
#remove(tSNETotal)

As with the correlation heatmap, we see that the samples from day 0 and 1 are clustered together while those from day 3 and 7 are clustered together.

t-SNE of Day 0 and Day 1

Now that we have seen that the samples from day 0 and 1 are closely related. Let’s try to take only those 2 in order to see if we can separate them from each other.

#tsne on Day 0/Day 1 cells only
indexDay01 <- which(data_day %in% c("Day 0","Day 1"))
group_one_data <- alive_tsne_data[indexDay01, ] #selecting only day 0 and day 1 cells
tsne_out <- Rtsne(dist(group_one_data)) # Run TSNE
# Show the objects in the 2D tsne representation
tsneDay01 <- qplot(x=tsne_out$Y[,1], y=tsne_out$Y[,2],color=data_day[indexDay01]) +
  labs(colour = "Cell Type", x = "tsne1", y = "tsne2") +
  ggtitle("t-SNE only day 0 and 1") +
  theme(text = element_text(size=14),plot.title = element_text(hjust = 0.5)) 
tsneDay01

#will not be needed: don't keep big data
remove(tsne_out)
#remove(tsneDay01)

The t-SNE can clearly not separate the 2 samples from each other. Finally let’s investigate whether they can be separated with a t-SNE in 3 dimension:

tsne_out <- Rtsne(dist(group_one_data), dims = 3) # Run TSNE in 3 dimensions
plot_ly(x=tsne_out$Y[,1],
        y=tsne_out$Y[,2],
        z=tsne_out$Y[,3],
        type="scatter3d",
        color=data_day[indexDay01],
        mode="markers") %>% 
  layout(title = 't-SNE day 0 and 1 3D')

#will not be needed: don't keep big data
remove(tsne_out)
remove(group_one_data)
#remove(tsneDay01)

We can still not distinguish both type of cells. This really shows how similar they are.

t-SNE of Day 0 and Day 1

Now let’s look at the cluster Day 3 and 7, to see if they can be distinguished.

#tsne on Day 0/Day 1 cells only
indexDay37 <- which(data_day %in% c("Day 3","Day 7"))
group_two_data <- alive_tsne_data[indexDay37, ] #selecting only day 0 and day 1 cells
tsne_out <- Rtsne(dist(group_two_data)) # Run TSNE
# Show the objects in the 2D tsne representation
tsneDay37 <- qplot(x=tsne_out$Y[,1], y=tsne_out$Y[,2],color=data_day[indexDay37]) +
  labs(colour = "Cell Type", x = "tsne1", y = "tsne2") +
  ggtitle("t-SNE only day 3 and 7") +
  theme(text = element_text(size=14),plot.title = element_text(hjust = 0.5)) 
tsneDay37

#will not be needed: don't keep big data
remove(tsne_out)
#remove(tsneDay37)

As before, the 2 type of cells cannot be distinguished with a t-SNE plot in 2D. They will probably not be distinguishable in 3D but let’s be sure:

tsne_out <- Rtsne(dist(group_two_data), dims = 3) # Run TSNE in 3 dimensions
plot_ly(x=tsne_out$Y[,1],
        y=tsne_out$Y[,2],
        z=tsne_out$Y[,3],
        type="scatter3d",
        color=data_day[indexDay37],
        mode="markers") %>% 
  layout(title = 't-SNE day 3 and 7 3D')

#will not be needed: don't keep big data
remove(tsne_out)
remove(group_two_data)

Not suprisingly, the second cluster cannot be seprated further with 3 dimensional t-SNE.

Fit linear model

Note that for fitting linear models, we cannot use the DATASEQ library as it requires counts (i.e. discrete values). We now have continuous values because of they way we normalized the data.

Voom

Mean-Variance trend

First let us visualize the mean-variance trend, to see. Note: here it is important that we use only the genes with no non zero values for every sample as we can see in paragraph 15.3 of (limma manual)[https://www.bioconductor.org/packages/devel/bioc/vignettes/limma/inst/doc/usersguide.pdf]. From the same article they say that as a rule of thumb we should keep only genes that have a

#working ... HOW TO USE DATADESEQ ON NORMALIZED => NOT COUNTS ANY MORE?? + IS IT WORTH IT ? cAN we round ?
#dataDESeq <- rlog(as.matrix(data), blind=FALSE )
#makes a convenient dge object for analysis
#t <- n_data[which(rowSums(n_data) > 300),] 
dge <- DGEList(counts=dataNZ, group = data_day)
# applies TMM normalization to dge
#isExpr <- rowSums(dge$counts) > 1000
#dge <- dge[isExpr,]
dge <- calcNormFactors(dge)
data_voomed <- voom(dge,design,plot=TRUE)

#will not be needed: don't keep big data
remove(dge)

The mean-variance trend is suprising: it doesn’t look flat enough. The first hyporthesis we had is that we there was a problem because we normalize by samples and then limma does it again. But not normailizing our data doesn’t change the curve of the plot. We then thought that maybe limma was doing something that required discrete data (as it is made for count data), our data isn’t discrete as we normailized it based on spike RNA. But rounding the data doesn’t change anything either. We will have to look deeper at the meaning of te curve of the plot and whether our data doesn’t follow the Voom assumptions. Note that the data may also be strange due to the nature of miRNA data, indeed in RNAseq we often have high value of counts but in miRNA the number of counts are much lower. This could cause the spike we see below x = zero (i.e. very low mean).

We also see that the mean variance trend seam to be decreasing very quickly and steadily, we can thus conclude that the experiment seem to be of low biological variation. For more information please see the “Removing heteroscedascity from count data” subpart of (this paper)[https://f1000research.com/articles/5-1408]

Here is an example of the data we have:

data_voomed[1:2,1:10]

An object of class "EList"
$targets
               group  lib.size norm.factors
HL60D0.C01_S1  Day 0 100862.30    0.9833482
HL60D0.C02_S2  Day 0  68768.90    0.9420345
HL60D0.C03_S3  Day 0  60552.35    1.0147531
HL60D0.C04_S4  Day 0 279868.19    0.9569073
HL60D0.C05_S5  Day 0  65361.77    1.1399012
HL60D0.C06_S6  Day 0  36823.33    0.8372045
HL60D0.C07_S7  Day 0 127944.17    0.8744118
HL60D0.C08_S8  Day 0  51630.86    0.9873263
HL60D0.C09_S10 Day 0 217824.43    1.0843357
HL60D0.C10_S11 Day 0  55444.19    1.1485957

$E
              HL60D0.C01_S1 HL60D0.C02_S2 HL60D0.C03_S3 HL60D0.C04_S4 HL60D0.C05_S5 HL60D0.C06_S6
hsa-let-7a-5p     13.374469     12.879195      12.91583     7.6516046      13.06635      13.48863
hsa-let-7a-3p      2.309527      2.862079      10.11151     0.8371754      10.14368      10.50555
              HL60D0.C07_S7 HL60D0.C08_S8 HL60D0.C09_S10 HL60D0.C10_S11
hsa-let-7a-5p     13.102009      13.73708       12.72418      13.310086
hsa-let-7a-3p      1.966402      10.01462        9.35818       9.870383

$weights
          [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]      [,8]      [,9]
[1,] 0.5885646 0.5004936 0.4749981 0.9426669 0.4901102 0.3899081 0.6533318 0.4453390 0.8348369
[2,] 0.2090121 0.1851257 0.1774882 0.2881312 0.1820762 0.1454568 0.2249386 0.1678079 0.2654087
         [,10]
[1,] 0.4582897
[2,] 0.1720934

$design
   HL60D0 HL60D1 HL60D3 HL60D7
1       1      0      0      0
2       1      0      0      0
3       1      0      0      0
4       1      0      0      0
5       1      0      0      0
6       1      0      0      0
7       1      0      0      0
8       1      0      0      0
9       1      0      0      0
10      1      0      0      0

Mean variance trend no log

Note: Seeing the strange plot above makes us think that maybe our count data doesn’t need to be log transformed. Let’s investigate the relationship between the mean and the variance of ou data.

#working ... HOW TO USE DATADESEQ ON NORMALIZED => NOT COUNTS ANY MORE?? + IS IT WORTH IT ? cAN we round ?
#dataDESeq <- rlog(as.matrix(data), blind=FALSE )
#takes only genes that are non zero in many samples
#t <- n_data[which(rowSums(n_data) > 300),] 
#tries plotting directly without voom
m <- rowMeans(dataNZ)
RowSigmaSqrt <- function(x) {
  (rowSums((x - rowMeans(x))^2)/(dim(x)[2] - 1))^(1/4)
}
v <- RowSigmaSqrt(dataNZ)
plotxy <- data.frame(mean=m,sigmaSqrt=v)
ggplot(plotxy, aes(x=mean, y=sigmaSqrt)) +
    geom_point(shape=1) +    # Use hollow circles
    geom_smooth(method=lm,   # Add linear regression line
                se=FALSE) + # Don't add shaded confidence region
  ylim(0, 7500) + 
  xlim(0,7500) + 
  ggtitle("Mean variance relationship not Voom") +
  theme(text = element_text(size=14),plot.title = element_text(hjust = 0.5)) 

From this plot we clearly see that our variance grows with fourth poqer of our mean. The reasoning behind this different relationship is still unclear to us (????????????? biologically), but it could indicate that we shouldn’t work with the default limma package and that we should maybe do the analysis manually (????????).

Fit linear model

fit <- lmFit(data_voomed, design)
fit <- contrasts.fit(fit, contrasts=contrastMatrix)
efit <- eBayes(fit)
topTable(efit, coef=colnames(coef(efit)))

plotSA(efit)

#will not be needed: don't keep big data
remove(fit)
remove(data_voomed)

Here there is a problem as we do not see the fitted line (?????).

Examining the number of DE genes

Let’s now look at differential expression levels of genes. Note that the adjusted p-value cutoff is 5%.

dt <- decideTests(efit)
summary(dt)

   HL60D1 - HL60D0 HL60D3 - HL60D0 HL60D7 - HL60D0 HL60D1 - HL60D3 HL60D1 - HL60D7
-1              28             325             498              87              41
0              705             376             206             292             221
1               11              43              40             365             482
   HL60D3 - HL60D7
-1              24
0              624
1               96

As we have already previously seen in the multiple plots of the previous section: most differential expressed genes are between day 0/1 and 3/7. We also note that a lot of genes in Day 3 and 7 are downregulated (compared to day 0 and 1).

Finding genes that are often differentially expressed compared to day 0

Let’s look at the genes that are always differentially expressed compared to day 0

# Finding genes that are not zero in every comparaison with day 0
commonDEgenesDay0 <- which(dt[,1]!=0 & dt[,2]!=0 & dt[,3]!=0)
names(commonDEgenesDay0)

 [1] "hsa-let-7a-3p"   "hsa-let-7b-3p"   "hsa-miR-1260a"   "hsa-miR-132-3p"  "hsa-miR-132-5p" 
 [6] "hsa-miR-143-3p"  "hsa-miR-17-3p"   "hsa-miR-181a-3p" "hsa-miR-181b-5p" "hsa-miR-181b-3p"
[11] "hsa-miR-199b-5p" "hsa-miR-19b-3p"  "hsa-miR-28-3p"   "hsa-miR-30d-5p"  "hsa-miR-5701"   
[16] "hsa-miR-7974"   

Let’s visualise it with a nice venn diagramm

vennDiagram(dt[,1:3], circle.col=c("turquoise", "salmon", "green"))

Finding genes that are often differentially between day 0/1 and day 3/7

Let’s look at the genes that are always differentially between day 0/1 and 3/7 as they seem to be the ones of interest.

# Finding genes that are not zero in every comparaison with day 0
commonDEgenesDay0137 <- which(dt[,2]!=0 & dt[,3]!=0 & dt[,4]!=0 & dt[,5]!=0)
names(commonDEgenesDay0137)

  [1] "hsa-let-7c-3p"     "hsa-let-7d-3p"     "hsa-let-7e-5p"     "hsa-let-7i-5p"    
  [5] "hsa-miR-101-3p"    "hsa-miR-106b-3p"   "hsa-miR-106b-5p"   "hsa-miR-1233-3p"  
  [9] "hsa-miR-1234-3p"   "hsa-miR-1236-3p"   "hsa-miR-125a-5p"   "hsa-miR-125b-2-3p"
 [13] "hsa-miR-1260a"     "hsa-miR-1268a"     "hsa-miR-1268b"     "hsa-miR-1271-5p"  
 [17] "hsa-miR-1271-3p"   "hsa-miR-1273d"     "hsa-miR-1273h-5p"  "hsa-miR-1275"     
 [21] "hsa-miR-1278"      "hsa-miR-1281"      "hsa-miR-1285-3p"   "hsa-miR-1287-5p"  
 [25] "hsa-miR-129-2-3p"  "hsa-miR-1291"      "hsa-miR-1303"      "hsa-miR-1304-3p"  
 [29] "hsa-miR-1306-3p"   "hsa-miR-1307-3p"   "hsa-miR-132-3p"    "hsa-miR-138-5p"   
 [33] "hsa-miR-143-5p"    "hsa-miR-143-3p"    "hsa-miR-145-5p"    "hsa-miR-146b-3p"  
 [37] "hsa-miR-149-5p"    "hsa-miR-150-5p"    "hsa-miR-151a-5p"   "hsa-miR-152-5p"   
 [41] "hsa-miR-1537-5p"   "hsa-miR-155-3p"    "hsa-miR-181a-3p"   "hsa-miR-181b-3p"  
 [45] "hsa-miR-181c-3p"   "hsa-miR-182-5p"    "hsa-miR-1827"      "hsa-miR-186-3p"   
 [49] "hsa-miR-190a-5p"   "hsa-miR-192-5p"    "hsa-miR-193a-5p"   "hsa-miR-193a-3p"  
 [53] "hsa-miR-195-5p"    "hsa-miR-196a-5p"   "hsa-miR-196b-5p"   "hsa-miR-197-5p"   
 [57] "hsa-miR-19b-3p"    "hsa-miR-200c-5p"   "hsa-miR-205-5p"    "hsa-miR-210-3p"   
 [61] "hsa-miR-2115-3p"   "hsa-miR-2116-3p"   "hsa-miR-212-3p"    "hsa-miR-218-5p"   
 [65] "hsa-miR-223-3p"    "hsa-miR-224-5p"    "hsa-miR-2277-3p"   "hsa-miR-2277-5p"  
 [69] "hsa-miR-2278"      "hsa-miR-23a-5p"    "hsa-miR-23a-3p"    "hsa-miR-25-3p"    
 [73] "hsa-miR-25-5p"     "hsa-miR-27a-5p"    "hsa-miR-29b-1-5p"  "hsa-miR-301b-5p"  
 [77] "hsa-miR-3074-3p"   "hsa-miR-3115"      "hsa-miR-3128"      "hsa-miR-3136-5p"  
 [81] "hsa-miR-3146"      "hsa-miR-3154"      "hsa-miR-3158-3p"   "hsa-miR-3173-3p"  
 [85] "hsa-miR-3177-3p"   "hsa-miR-3180-5p"   "hsa-miR-3181"      "hsa-miR-3198"     
 [89] "hsa-miR-328-3p"    "hsa-miR-338-5p"    "hsa-miR-33b-3p"    "hsa-miR-34a-5p"   
 [93] "hsa-miR-3529-3p"   "hsa-miR-3607-5p"   "hsa-miR-3607-3p"   "hsa-miR-3614-3p"  
 [97] "hsa-miR-3648"      "hsa-miR-365a-3p"   "hsa-miR-365b-3p"   "hsa-miR-3662"     
[101] "hsa-miR-3677-3p"   "hsa-miR-3680-5p"   "hsa-miR-3680-3p"   "hsa-miR-3691-5p"  
[105] "hsa-miR-378b"      "hsa-miR-378g"      "hsa-miR-3912-5p"   "hsa-miR-3913-5p"  
[109] "hsa-miR-3917"      "hsa-miR-3922-3p"   "hsa-miR-3940-3p"   "hsa-miR-3940-5p"  
[113] "hsa-miR-3942-5p"   "hsa-miR-3944-3p"   "hsa-miR-409-3p"    "hsa-miR-422a"     
[117] "hsa-miR-424-3p"    "hsa-miR-425-5p"    "hsa-miR-4284"      "hsa-miR-4421"     
[121] "hsa-miR-4426"      "hsa-miR-4429"      "hsa-miR-4434"      "hsa-miR-4435"     
[125] "hsa-miR-4449"      "hsa-miR-4454"      "hsa-miR-4477b"     "hsa-miR-4484"     
[129] "hsa-miR-450a-2-3p" "hsa-miR-4510"      "hsa-miR-4515"      "hsa-miR-4518"     
[133] "hsa-miR-4523"      "hsa-miR-455-5p"    "hsa-miR-4660"      "hsa-miR-4662a-5p" 
[137] "hsa-miR-4664-5p"   "hsa-miR-4685-3p"   "hsa-miR-4707-5p"   "hsa-miR-4723-5p"  
[141] "hsa-miR-4743-3p"   "hsa-miR-4744"      "hsa-miR-4745-5p"   "hsa-miR-4745-3p"  
[145] "hsa-miR-4747-3p"   "hsa-miR-4766-3p"   "hsa-miR-4772-3p"   "hsa-miR-4777-3p"  
[149] "hsa-miR-4779"      "hsa-miR-4782-3p"   "hsa-miR-4785"      "hsa-miR-4791"     
[153] "hsa-miR-4796-5p"   "hsa-miR-4797-3p"   "hsa-miR-4802-3p"   "hsa-miR-486-5p"   
[157] "hsa-miR-495-3p"    "hsa-miR-497-5p"    "hsa-miR-5001-5p"   "hsa-miR-5004-3p"  
[161] "hsa-miR-5094"      "hsa-miR-5188"      "hsa-miR-548a-3p"   "hsa-miR-548aa"    
[165] "hsa-miR-548ag"     "hsa-miR-548am-5p"  "hsa-miR-548ar-3p"  "hsa-miR-548au-3p" 
[169] "hsa-miR-548aw"     "hsa-miR-548ay-5p"  "hsa-miR-548b-5p"   "hsa-miR-548d-5p"  
[173] "hsa-miR-548f-3p"   "hsa-miR-548n"      "hsa-miR-548s"      "hsa-miR-548t-5p"  
[177] "hsa-miR-548t-3p"   "hsa-miR-550b-2-5p" "hsa-miR-552-5p"    "hsa-miR-5585-5p"  
[181] "hsa-miR-5586-5p"   "hsa-miR-5684"      "hsa-miR-5690"      "hsa-miR-570-3p"   
[185] "hsa-miR-573"       "hsa-miR-581"       "hsa-miR-582-5p"    "hsa-miR-597-3p"   
[189] "hsa-miR-6134"      "hsa-miR-617"       "hsa-miR-628-3p"    "hsa-miR-629-5p"   
[193] "hsa-miR-639"       "hsa-miR-6500-3p"   "hsa-miR-6509-5p"   "hsa-miR-6513-3p"  
[197] "hsa-miR-663a"      "hsa-miR-663b"      "hsa-miR-6724-5p"   "hsa-miR-6726-3p"  
[201] "hsa-miR-6727-5p"   "hsa-miR-6729-3p"   "hsa-miR-6730-5p"   "hsa-miR-6733-3p"  
[205] "hsa-miR-6734-3p"   "hsa-miR-6734-5p"   "hsa-miR-6735-3p"   "hsa-miR-6737-3p"  
[209] "hsa-miR-6749-3p"   "hsa-miR-6753-5p"   "hsa-miR-6754-3p"   "hsa-miR-6762-3p"  
[213] "hsa-miR-6763-5p"   "hsa-miR-6764-3p"   "hsa-miR-6764-5p"   "hsa-miR-6774-3p"  
[217] "hsa-miR-6777-3p"   "hsa-miR-6781-3p"   "hsa-miR-6781-5p"   "hsa-miR-6787-3p"  
[221] "hsa-miR-6790-3p"   "hsa-miR-6793-3p"   "hsa-miR-6797-3p"   "hsa-miR-6798-3p"  
[225] "hsa-miR-6803-3p"   "hsa-miR-6812-3p"   "hsa-miR-6813-5p"   "hsa-miR-6818-3p"  
[229] "hsa-miR-6819-3p"   "hsa-miR-6833-3p"   "hsa-miR-6837-5p"   "hsa-miR-6844"     
[233] "hsa-miR-6849-3p"   "hsa-miR-6852-5p"   "hsa-miR-6857-3p"   "hsa-miR-6858-5p"  
[237] "hsa-miR-6859-3p"   "hsa-miR-6862-3p"   "hsa-miR-6868-3p"   "hsa-miR-6873-3p"  
[241] "hsa-miR-6885-3p"   "hsa-miR-6889-3p"   "hsa-miR-6891-5p"   "hsa-miR-6891-3p"  
[245] "hsa-miR-6894-3p"   "hsa-miR-7111-3p"   "hsa-miR-7112-3p"   "hsa-miR-7114-3p"  
[249] "hsa-miR-7641"      "hsa-miR-766-5p"    "hsa-miR-7705"      "hsa-miR-7854-3p"  
[253] "hsa-miR-8485"      "hsa-miR-874-3p"    "hsa-miR-877-5p"    "hsa-miR-885-5p"   
[257] "hsa-miR-9-3p"      "hsa-miR-92a-1-5p"  "hsa-miR-92a-3p"    "hsa-miR-937-5p"   
[261] "hsa-miR-940"       "hsa-miR-942-3p"    "hsa-miR-99a-3p"   

Let’s visualise it with a nice venn diagramm

vennDiagram(dt[,2:5], circle.col=c("turquoise", "salmon", "green","yellow"))